clc;clear;
% syms t
% q=1000;
% pM=101.325*10^3;
% pv=2.318*10^3;
% sigma=0.0728;
% u=1002*10^(-6);
% ye=1.4;
% A=0.4;
% R0=0.26*10^(-3);
% y=R0;
% Dy=130.18;
% D=(-3/2*Dy^2-1/q*(pM-pv-(A*pM-pv+2*sigma/R0)*(R0/y)^(3*ye))-(4*u/q)*(Dy/y)-2*sigma/(q*y))/y;
% y0=dsolve('D2y*y+3/2*Dy^2+1/q*(pM-pv-(A*pM-pv+2*sigma/R0)*(R0/y)^(3*ye))+(4*u/q)*(Dy/y)+2*sigma/(q*y)','y(0)=R0,Dy(0)=0','t')  %用dsolve求解，表达式，初值，指定自变量
% t=0:0.05*(-3):1.5*(-3); %求解范围划定
% y0=eval(y0);             %得到数值解
% plot(t,y0)

%%dsolve求解不了
R0=0.26*10^(-3);
[t,R]=ode15s(@fx,[0,4.5*10^(-3)],[R0,130.18]);  %第一项需要定义函数里面包含函数表达式，第二项是t的范围，第三项是初值
plot(t,R(:,1))                                  %画图
xlabel("t(s)")
ylabel("R(m)")

%求加速度，因为前面的R只能输出半径和速度，这里重新计算
q=1000;
pM=101.325*10^3;
pv=2.318*10^3;
sigma=0.0728;
u=1002*10^(-6);
ye=1.4;
A=0.4;
R0=0.26*10^(-3);
%——————————参数值带入
S=size(R,1);      %确定R的行数
R1=R(:,1);        %半径
R2=R(:,2);        %速度
res=zeros(S,1);   %加速度保存在res里
for i=1:S
res(i)=(-3/2*R2(i)^2-1/q*(pM-pv-(A*pM-pv+2*sigma/R0)*(R0/R1(i))^(3*ye))-(4*u/q)*(R2(i)/R1(i))-2*sigma/(q*R1(i)))/R1(i);
end
